Michael Ruggiero
# Function imports
import numpy as np
import pandas as pd
import networkx as nx
import osmnx as ox
import geopandas as gpd
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import matplotlib.colors as colors
ox.config(use_cache=True, log_console=True)
ox.__version__
%config InlineBackend.figure_format = 'retina'
%matplotlib inline
1) Build a NLP model identifies road closures based on in social media
2) Use Google maps to chart travel based on different times of day
3) Translate "bad traffic areas" as flaged map locations
4) Simulate a natural disaster that destroys roads
5) Reflect the affected roads into a map in red color
6) Produce a map of valid rescue roads/escape routes
7) Optimize dispatch for police and rescue
# get a graph for some city
place = {'city' : 'Medford',
'state' : 'MA',
'country' : 'USA'}
medford = ox.graph_from_place(place, network_type='bike')
fig, ax = ox.plot_graph(medford)
# what sized area does our network cover in square meters?
medford_proj = ox.project_graph(medford)
nodes_proj = ox.graph_to_gdfs(medford_proj, edges=False)
graph_area_m = nodes_proj.unary_union.convex_hull.area
print(graph_area_m)
# show basic stats about the network
medford_stats = ox.basic_stats(medford_proj,
area=graph_area_m,
clean_intersects=True,
circuity_dist='euclidean')
# medford_stats = ox.extended_stats(medford_proj, ecc=True, bc=True, cc=True)
pd.DataFrame(medford_stats).set_index("streets_per_node_counts")
# Edge and Node projection
nodes_med, edges_med = ox.graph_to_gdfs(medford, nodes=True, edges=True)
nodes_med.head(2)
nodes_med.shape
#Functions built for this project
# %load graph_functions.py
import graph_functions as gf
medford_data = gf.node_roader(medford)
medford_data.keys()
edges_med, nodes_med, medford = medford_data["edges"], medford_data["nodes"], medford_data["graph"]
nodes_med.major_inter.value_counts()
#Long/Lat information
nodes_med[["x","y"]].describe()
nodes_med.head(2)
#A function to generate a random lat/long point
gf.random_point()
random_point_A = gf.random_point()
random_point_B = gf.random_point()
# get the nearest network node to each random point
orig_node = ox.get_nearest_node(medford, random_point_A)
dest_node = ox.get_nearest_node(medford, random_point_B)
# find the route between these nodes then plot it
route = nx.shortest_path(medford, orig_node, dest_node, weight='length')
fig, ax = ox.plot_graph_route(medford, route, node_size=0)
# Length of route in meters
nx.shortest_path_length(medford, orig_node, dest_node, weight='length')
# Absolute Distance between two nodes (as the crow flies)
ox.great_circle_vec(medford.node[orig_node]['y'], medford.node[orig_node]['x'],
medford.node[dest_node]['y'], medford.node[dest_node]['x'])
d_location = nodes_med[["y","x"]].sample(1).values[0]
print(d_location)
#returns a graphs object of a disater
disaster = gf.disaster_generator(nodes_med, location_point = d_location, radius = 2500)
nodes_dis, edges_dis = ox.graph_to_gdfs(disaster)
#Set Colors
ec = ['lightcoral' if i in disaster.edges() else 'green' for i in medford.edges()]
nc = ['red' if i in disaster.nodes() else 'blue' for i in medford.nodes()]
#Plot energency grid
ox.plot_graph(medford, node_size=15, node_color = nc ,edge_color=ec)
nx.degree_histogram(disaster)
medford_disaster = gf.road_kill(disaster,
nodes_med,
edges_med,
kill_percentage = .8 )
# get the boundary polygons for neighboring cities, save as shapefile, project to UTM, and plot
place_names = ['Medford, MA, USA',
'Arlington, MA, USA',
'Somerville, MA, USA',
'Malden, MA, USA',
# 'Melrose, MA, USA',
'Winchester, MA, USA',
'Stoneham, MA, USA',
'Everett, MA, USA',
]
medford_area = ox.gdf_from_places(place_names)
fig, ax = ox.plot_shape(medford_area)
# highlight one-way roads
oneway = ox.graph_from_place(place, network_type='drive')
ec = ['r' if data['oneway'] else 'b' for u, v, key, data in oneway.edges(keys=True, data=True)]
fig, ax = ox.plot_graph(oneway, node_size=0, edge_color=ec, edge_linewidth=1.5, edge_alpha=0.5)
#Neighboor hoods buffered
neigboors_buffered = ox.gdf_from_places(place_names, gdf_name='neighboors', buffer_dist=250)
fig, ax = ox.plot_shape(neigboors_buffered, alpha=0.7)
# DF builder of start/end nodes with long/lat and times
# Already done, commenting out to avoid clutter in notebook
router = gf.random_zone_picker(nodes_med, 10)
#router.to_csv("more_zones.csv")
router.head(3)
# network from address, including only nodes within 5.5km along the network from city hall
neighbor = ox.graph_from_address(address='85 George P. Hassett Drive,\
Medford, MA 02155',
distance=5500,
distance_type='network',
network_type='walk')
# you can project the network to UTM (zone calculated automatically)
neighbor_projected = ox.project_graph(neighbor)
# Edge and Node projection
nodes_area, edges_area = ox.graph_to_gdfs(neighbor_projected,
nodes=True,
edges=True)
nc = ['blue' if i in medford.nodes() else 'green' for i in neighbor.nodes()]
ox.plot_graph(neighbor,
node_size=15,
node_color = nc);
neighbor_data = gf.node_roader(neighbor)
black = neighbor_data["nodes"]
black[black.color == "black"].head(2)
neighbor_data.keys()
nodes_area, edges_area, neighborhood = neighbor_data["nodes"], neighbor_data["edges"], neighbor_data["graph"]
for i in neighbor_data['major_intersections']:
print("Number of major nodes: ",
len(neighbor_data['major_intersections'][i]),
"\t with degree ",
i,)
neighbor_data['major_intersections'].keys()
# find the route between these nodes then plot it
route = nx.shortest_path(neighbor, 3809869827, 65339393, weight='length')
fig, ax = ox.plot_graph_route(neighbor, route, node_size=0)
_ = gf.city_node(nodes_med,"medford",nodes_area)
nodes_area.head()
gmap, complete_node, node_occurances = gf.directionsToDataframe("10.Google_directions_api_all6k.csv")
gmap.head()
google_top = gf.traffic_node(nodes_area, node_occurances)
g_list = list(google_top[google_top.traffic_importance > 887].osmid.values)
google_top[google_top.traffic_importance > 887]
# Edge and Node projection
nodes_area, edges_area = ox.graph_to_gdfs(neighbor_projected,
nodes=True,
edges=True)
g_size = [25 if i in g_list else 0 for i in neighbor.nodes()]
ox.plot_graph(neighbor,
node_size=g_size,
node_color = "red", edge_alpha= .5);
import time
import pickle
from scipy.optimize import linear_sum_assignment
#Populate police station
"""
Medford Police:
100 Main St, Medford, MA 02155
42.415811 -71.110152
"""
police_node = gf.nodefinder(nodes_med,
lat = 42.415811,
lng = -71.110152)
def monte_simulator(p_e_samples, police_location, simulations, name = ""):
main_start = time.time()
police_node = gf.nodefinder(nodes_med,
lat = police_location[0],
lng = police_location[1])
police_station = [police_node] * p_e_samples
p_e_samples *= 2
monte_simulation = {}
#Iteration though simulation
for sim in range(0,simulations):
#Initialization
start = time.time()
print("Starting simulation: {} at time: {}".format(sim,round(start - main_start,2)))
monte_dis = {}
#Generates disaster and radius
disaster_location = nodes_med.sample(1)
monte_dis["dis_location"] = disaster_location
monte_dis["dis_rad"] = np.random.uniform(low = .4, high = .9) * 4000
monte_dis["disaster"] = gf.disaster_generator(nodes_med,
location_point = disaster_location[["y","x"]].values[0],
radius = monte_dis["dis_rad"],
plotter = 0)
monte_dis["nodes_effected"] = len(monte_dis["disaster"].nodes())
monte_dis["roads_effected"] = len(monte_dis["disaster"].edges())
#Randomly destroys roads in disaster radius
monte_dis["remain_per"] = np.random.uniform(low = .35, high = 1)
monte_dis["remaining"] = gf.road_kill(monte_dis["disaster"],
nodes_area, #These come from the neighborhood
edges_area,
kill_percentage = monte_dis["remain_per"],
plotter = 0)
#Make dictionaries for emergencies and police
e_list = np.random.choice(list(monte_dis["disaster"].nodes()),
replace=False,
size = p_e_samples)
p_list = list(nodes_med.sample(int(p_e_samples / 2)).osmid.values)
p_list.extend(police_station)
monte_dis["emergency"] = {}
monte_dis["patrol"] = {}
for i in range(0,p_e_samples):
monte_dis["patrol"]["officer_{}".format(i)] = p_list[i]
monte_dis["emergency"]["emergency_{}".format(i)] = e_list[i]
rows = list(monte_dis["patrol"].keys())
columns = list(monte_dis["emergency"].keys())
#Distance Matrix
d_matrix = pd.DataFrame(columns = columns,
index = rows)
for i in range(len(rows)):
for j in range(len(columns)):
try:
path = nx.shortest_path_length(monte_dis["remaining"],
monte_dis["patrol"][rows[i]],
monte_dis["emergency"][columns[j]],
weight='length')
except:
path = 10_000
d_matrix.loc[rows[i], columns[j]] = path
monte_dis["distance_matrix"] = d_matrix.astype(float)
#Use Hungarian algorithm for linear sum assignment
lsa = linear_sum_assignment(monte_dis["distance_matrix"].to_numpy())
adj_matrix = np.matrix(np.zeros((10, 10)))
#optimized distances
od_list = []
for i in range(len(lsa[0])):
od = monte_dis["distance_matrix"].iloc[lsa[0][i],lsa[1][i]]
adj_matrix[lsa[0][i],lsa[1][i]] = od
od_list.append(od)
dispatch = {"officer_{}".format(lsa[0][i]):["emergency_{}".format(lsa[1][i]),
od_list[i]] for i in range(len(lsa[0]))}
monte_dis["ad_matrix"] = adj_matrix
monte_dis["dispatch"] = dispatch
monte_dis["dispatch_distance"] = sum(od_list)
monte_simulation[sim] = monte_dis
if sim % 200 == 0:
print("pickle dump")
pickle.dump(monte_simulation, open( "monte_simulation_{}_{}.p".format(name,sim), "wb" ))
end = time.time()
print("Ending simulation: {}. Total {} seconds for run\n".format(sim,round(end - start,2)))
return monte_simulation
monte_simulation = monte_simulator(5,[42.415811,-71.110152], 5, "demo")
monte_simulation.keys()
monte_simulation[0].keys()
monte_simulation[0]['dispatch']
monte_simulation[2]['dispatch']
monte_simulation[2]["ad_matrix"]
size_list = range(len(monte_simulation))
df = pd.concat([monte_simulation[i]["dis_location"].copy() for i in size_list]).reset_index(drop=True)
extra_data = ['dis_rad', 'nodes_effected', 'roads_effected', 'remain_per', 'dispatch_distance']
for i in range(len(monte_simulation)):
for j in extra_data:
df.loc[i, j] = monte_simulation[i][j]
df.columns
X = df[['osmid','x', 'y','dis_rad', 'nodes_effected', 'roads_effected',
'remain_per', 'dispatch_distance']]
X.head()